An integer GARCH model for a Poisson process with time-varying zero-inflation

A serially dependent Poisson process with time-varying zero-inflation is proposed. Such formulations have the potential to model count data time series arising from phenomena such as infectious diseases that ebb and flow over time. The model assumes that the intensity of the Poisson process evolves according to a generalized autoregressive conditional heteroscedastic (GARCH) formulation and allows the zero-inflation parameter to vary over time and be governed by a deterministic function or by an exogenous variable. Both the expectation maximization (EM) and the maximum likelihood estimation (MLE) approaches are presented as possible estimation methods. A simulation study shows that both parameter estimation methods provide good estimates. Applications to two real-life data sets on infant deaths due to influenza show that the proposed integer-valued GARCH (INGARCH) model provides a better fit in general than existing zero-inflated INGARCH models. We also extended a non-linear INGARCH model to include zero-inflation and an exogenous input. This extended model performed as well as our proposed model with respect to some criteria, but not with respect to all.


Introduction
The standard Poisson point process, which assumes statistical independence between observations, is not suitable for modelling time series of counts that display serial dependence. For example, counts of infectious disease occurrences can be considered as dependent on previous counts because of the infectious nature of the illness. One way to address this deficiency is to define a time series where the count at a given time is generated by a Poisson distribution whose intensity parameter is dependent on past counts and past intensities. Observe that for a discrete time process defined in such a manner with a linear dependent structure and observed at equally spaced time points, the intensity parameter at a given time is the mean count at that time conditional on the past information. Since many count data processes have an excess of zeros than what is possible by the underlying probability distribution, that is they exhibit zeroinflation, modelling serial dependence is not sufficient. Existing models for count data that accommodate serial dependence and zero-inflation, however, do not allow for the zero-inflation probability to vary over time cyclically or be driven by an exogenous variable. This a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 becomes a handicap when modelling count data time series, such as infectious disease or death counts that vary seasonally or with respect to an external factor that varies with time. To address this shortcoming, we propose the Time-Varying Zero-Inflated Poisson integer-valued GARCH model (TVZIP-INGARCH), which is based on a zero-inflated Poisson process whose intensity parameter satisfies the integer-valued GARCH (INGARCH) formulation introduced by Ferland, Latour, and Oraichi [1]. It is also a generalization of the model introduced by Zhu [2] which assumes constant zero-inflation. The proposed model can accommodate cases where the zero-inflation probability is driven by a deterministic function of time, such as a sinusoidal wave, or by exogenous variables. Before we introduce this model formerly, we will discuss some relevant INGARCH type count data time series models that have been proposed in the past. Note that we extended an existing non-linear INGARCH model without zero-inflation to include a zero-inflation component as well as an exogenous input variable. The performance of this new formulation and several existing INGARCH type models with constant zero-inflation were compared to that of the proposed model with respect to model fit on two real-life data sets.
Rydberg and Shephard [3] were the first to propose a count data time series model that account for serial dependence. In their model, the current conditional mean is a linear function of both the observed count and the conditional mean at the pervious time point. Similar models have been proposed by several other authors and these are discussed in Chapter 4 of the book by Kedem and Fokianos [4]. Heinen [5] generalized the lag one model of Rydberg and Shephard to include an arbitrary number of lags for both the past counts and past conditional means and named it the Autoregressive Conditional Poisson model with lags p and q (ACP (p, q)). The formulation of this model resembles that of the generalized autoregressive conditional heteroscedastic (GARCH) model of Bollerslev [6], but models the conditional mean rather than the conditional variance. Heinen derived the properties of his model only for the ACP (1, 1) case. The general case was investigated by Ghahramani and Thavaneswaran [7] who referred to the Heinen paper as the origin of the ACP model. Independently, Ferland et al. [1] proposed what was termed the Integer-valued GARCH (INGARCH) process, which is essentially the same as the ACP model of Heinen. The INGARCH model of order (p, q) is defined as follows: where fX t : t 2 Zg is the count process, λ t defines the conditional mean of X t given the sigmafield F tÀ 1 generated by X l and λ l values for l<t. The model is defined with the parameter constraints a 0 > 0; a i � 0; b j � 0; i ¼ 1; . . . ; p; j ¼ 1; . . . ; q; p � 1; q � 0; for positivity of λ t and for stationarity of X t . In addition to presenting the model, Heinen [5] also derived the stationarity conditions, covariance functions, and addressed the problem of maximum likelihood estimation (MLE) of the parameters. As mentioned earlier, Zhu [2] extended this model to accommodate zero-inflation. In the following literature review, we will focus on additional work done on INGARCH models without zero-inflation first before moving onto the zero-inflated models. Weiß [8] extended the previous results on the class of INGARCH models and derived a set of Yule-Walker type equations for the autocorrelation function for the general INGARCH case. Several important theoretical contributions for both the linear and non-linear Poisson autoregressive models were made by Fokianos et al. [9]. For the model where the conditional mean is a linear function of the past conditional mean and past count, they proved that the maximum likelihood estimates are asymptotically normal. Note that the model they considered is INGARCH (1, 1), but the authors refrain from calling it as such and labelled it as Poisson autoregression because of their stated belief that the GARCH moniker should be reserved for formulations that model variance. In the non-linear case where the conditional mean is a non-linear function of its past values and past counts, the authors established that the intensities of the Poisson distribution at each time point form a geometrically ergodic Markov chain under some general assumptions. Fokianos and Tjøstheim [10] further extended these results and showed that in the log-linear Poisson autoregression case, the maximum likelihood estimates are asymptotically normal, and that covariance matrix of the parameter estimates are consistent. Chen and Lee [11] extended the log-linear Poisson INGARCH formulation by introducing linear effects exogenous variables to the equation for the logarithm of the conditional mean and proposed the log-linear Poisson INGARCHX and negative binomial INGARCHX models. Here X denotes a model with exogenous inputs. The authors suggest that the forecasting performance of the model can be improved by the addition of covariates to the model, and that it helps to properly analyse the causal relationships between exogenous factors and the count time series. The softplus INGARCH model, a non-linear version of INGARCH, was proposed by Weiß, Zhu, and Hoshiyar [12]. This model utilizes the softplus link function as an alternative to the log-link function and thus allowing for negative autocorrelations. We extended this model to include zero-inflation and accommodate an exogenous variable as input. Details of this extension are given in a later section.
Negative binomial (NB), Generalized Poisson (GP), and Double Poisson (DP) are well known discrete distributions that can also be used as an alternative to the Poisson process Zhu [13]. A negative binomial INGARCH model (NB-INGARCH), which is an alternative to the Poisson INGARCH formulation, was proposed and the stationary conditions and the autocorrelation function of the process were obtained by Zhu [14]. Ye, Garcia, Pourahmadi, and Lord [15] allowed the negative binomial INGRACH model from Zhu [14] to incorporate covariates, so that the relationship between a time series of counts and correlated external factors could be properly modelled. Xu, Xie, Goh, and Fu [16] adapted the INARCH formulation and proposed a new dispersed INGARCH (DINARCH) model to handle the conditional overdispersion and underdispersion. Xu et al. [16] also mention that when the dispersion parameter is not constant this model coincides with that of Zhu [14] without the moving average terms.
The idea of modelling the zero-inflated probability in a count data process as a function of covariates can be dated back to Lambert [17]. Zhu [2] while introducing the zero-inflated Poisson, also presented a zero-inflated negative binomial integer-valued GARCH model, but with the assumption of constant zero-inflation probability. In this paper two types of negative binomial distributions are discussed, namely Negative Binomial 1 (NB1), and Negative Binomial 2 (NB2). Chen and Lee [18] investigated the zero-inflated generalized Poisson autoregressive (ZIGP-AR) model of Lee, Lee, and Chen [19] and proposed the zero-inflated generalized Poisson INGARCH (ZIGP-INGARCH) which allows for structural breaks. In the ZIGP-IN-GARCH model zero-inflation is introduced to the generalized Poisson INGARCH (GP-INGARCH) model of Zhu [20]. Gonçalves, Lopes, and Silva [21] introduced a new class of zero-inflated INGARCH models that included general compound Poisson deviates and named it as the zero-inflated compound Poisson INGARCH (ZICP-INGARCH) process. This model can include both zero-inflated Poisson and zero-inflated negative binomial INGARCH models of Zhu [2].
More recent developments in the INGARCH literature are as follows. Xiong and Zhu [22] and Li, Chen, and Zhu [23] considered the robust estimation methods for INGARCH models.
Liu, Zhu, and Zhu [24] generalized the range of observations from infinite to categorical. Cui, Li and Zhu [25] and Xu and Zhu [26] generalized the INGARCH models from non-negative integer-valued to the integer-valued cases. Lee, Kim and Seok [27] introduced one parameter exponential family INGARCH model with zero-inflation and it was named as ZIEF-IN-GARCH. This new model can accommodate above mention zero-inflated versions of INGARCH models. In addition to that they developed residuals based CUSUM tests which can be used to examine the change points for the proposed model. A summary of various count time series published recently can be found in Davis et al. [28].
None of the above models allow for a time-varying zero-inflation probability. It is imperative that such accommodation should be made because empirical count data time with large number of zero counts can display strong cyclical behavior or seasonality with respect to the observed zero values. Ignoring this time-varying property of the zero-inflation parameter decreases the predictive performance of the model. Recognizing this, Yang [29], discussed the importance of modelling zero-inflation probability as a time-varying function. In the above article, it was assumed that both the zero-inflation and the intensity parameter are driven by the linear combination of past observations of exogenous variables and connects them to the conditional mean of the count via a log-link function. A recently introduced approach to modelling time-varying zero-inflation and the intensity parameter is the adaptive log-linear zero-inflated generalized Poisson model proposed by Xu et al. [30]. The authors assumed that the counts are from a zero-inflated generalized Poisson distribution, with the logarithm of the intensity propagated through a GARCH type model augmented with additional terms of exogenous variables and associated coefficients. All model parameters are assumed to be timedependent. If this time dependence and the exogenous variable are removed, then the model becomes the log-linear Poisson autoregressive formulation proposed by Fokianos and Tjøstheim [10]. If the time dependence and zero-inflation is removed, then it becomes the loglinear INGARCHX model of Chen and Lee [11]. In applying to monthly crime data from New South Wales, Australia, the authors assumed constant coefficients over a local interval at any given time point t (with the interval adaptively selected from among a prechosen set of nested intervals) and estimated the parameters separately using data from each such interval. In their approach, the shifting of the intervals allows for parameters to change from one interval to another and thus allowing for the zero-inflation as well as other model parameters to vary with time. Note that a version of this model with a time varying exogenous input but non-time varying zero-inflation was included in the comparison of our proposed model against several constant zero-inflation models in an empirical setting.
We propose a different approach based on a generalization of the model proposed by Zhu [2]. In our formulation, it is the zero-inflation probability, rather than the intensity of the Poisson process, that is allowed to be governed by exogenous variables. We also assume that the INGARCH model parameters remain constant over time, allowing for a more parsimonious model that is relatively easier to estimate. While the model proposed by Xu et al. [30] is more flexible, the approach we propose is presented as a simpler and practical alternative, albeit a less sophisticated one. As mentioned earlier, the proposed model can also accommodate the case where the zero-inflation probability is driven by a deterministic function of time or driven by an exogenous variable. In addition, the intensity parameter of the Poisson process is assumed to vary dynamically through a GARCH type model. Thus, the INGARCH part of the proposed model can be viewed as observation driven, in the sense that recursive substitutions can be employed to show that the current intensity of the process conditional on the past is a linear function of past observations and past intensities. Also note that we model the conditional mean rather than the logarithm of the conditional mean as is done in Xu et al. [30].
The remainder of this paper is organized as follows: the first, the Time-Varying Zero-Inflated Poisson INGARCH model (TVZIP-INGARCH) is introduced. Two cases, namely a deterministic cyclically varying zero-inflation component and a model with the zero-inflation parameter driven by an exogenous set of stochastic variables are discussed. A section on parameter estimation procedures is presented next, which is followed by the results of a simulation study. Following which is a section that introduces extensions to the softplus INGARCH model. The performance of this extended model together with additional zero-inflated INGARCH type models are compared to the proposed model in the next section. This section also presents the results of fitting the above models to two empirical data sets. A discussion and conclusions are presented in the final section of this paper.

The time-varying zero-inflated INGARCH model
As Zhu [2] stated, the probability mass function (pmf) of a zero-inflated Poisson model with parameter vector (λ,ω), with X representing the count, can be written in the following form: Zhu [2] presented the mean and the variance of the distribution as follows: EðXÞ ¼ lð1 À oÞ and VarðXÞ ¼ lð1 À oÞð1 þ loÞ > EðXÞ for 0 < o < 1: Moving on to define the proposed time-varying zero-inflated INGARCH model, assume that fX t : t 2 Zg is a discrete time series of count data, and F tÀ 1 is the sigma field generated by fX l ; l l : l � t À 1g. The conditional distribution of X t given F tÀ 1 and ω t , is assumed to be a zero-inflated Poisson (ZIP) with parameter vector (μ t ,ω t ) where m t ¼ ð1 À o t Þl t . Then, The dynamic propagation of the conditional mean (μ t ) of the zero-inflated Poisson process is defined by: with the intensity parameter λ t of the Poisson process formulated as: where a 0 > 0; a i � 0; b j � 0; i ¼ 1; 2; 3; . . . p; j ¼ 1; 2; 3; . . . ; q; p � 1; q � 1; and t 2 Z: Furthermore, o t ¼ gðV t ; ΓÞ 2 ð0; 1Þ8t 2 Z; is a function of variables, propagating over time, which is used to model the time-varying zero-inflation. Note that elements of the vector V t may consist of stochastic exogenous variables that vary with time, or it may be a scaler equal to time t. In addition, Γ denotes a vector of parameters. It is assumed that 0<ω t <1 for all t 2 Z: Note that Fokianos et al. [9] defined their linear Poisson autoregressive model for t 2 N instead of t 2 Z; with constant initial conditions. The model in (2) can also be defined in a similar manner. It is this alternative formulation that was used in our simulation study.
The above model is denoted by TVZIP-INGARCH (p, q). If p>0 and q = 0, then the model becomes a TVZIP-INARCH model with order p, denoted by TVZIP-INARCH (p). The conditional mean and conditional variance of X t given F tÀ 1 and ω t are specified by the following equations: The conditional variance to conditional mean ratio, or the dispersion ratio, of TVZIP distribution is: Note that results (3) and (4) can be derived using the definition of the conditional mean and variance and applying standard procedures utilized in deriving similar quantities related to GARCH type models.
The result in (4) indicates that TVZIP-INGARCH (p, q) can be used to model integer-valued time series with overdispersion if the values of ω t and λ t are uniformly bounded below by positive constants.
Case 1: Zero-inflation driven by a deterministic function of time. In Case 1, it is assumed that the zero-inflation function ω t =g(V t ,Γ) is such V t is a scaler equal to t. For example, we may assume the function g to be defined as follows: where s is the seasonal length, and Γ ¼ fðA; B; CÞ T jA; B; C 2 Rg.
As mentioned above, the time-varying zero-inflation function ω t = g(V t ,Γ) should always be bounded between zero and one. The range of values for A,B, and C in (5) that are needed to satisfy the above criterion are derived in S1 Appendix. Note that herein a simple example is used, where the function g consists of a sine function and a cosine function of equal period, but g could be any other function of time that, with proper selection of parameters, can be bounded between zero and one.
Case 2: Zero-inflated function driven by an exogenous variable. The proposed model also accommodates the case where the zero-inflation probability is determined by one or more exogenous variables. In this case g(V t ,Γ) is considered a function, with the interval (0, 1) as its range, of the vector of exogenous variables V t One example of g is the logistic function. Note that V t can be a scaler seasonal autoregressive time series, a vector seasonal time series, or a scaler or vector time series that varies non-seasonally. For illustrative purposes, consider the case where V t is a scaler purely seasonal autoregressive time series, denoted by V t , with period s, g the logistic function, and ε t is a white noise error term. Then we can write, with

Estimation procedure
The use of both the Expectation Maximization (EM) algorithm and Maximum Likelihood (ML) method to estimate the model parameters were developed for the general TVZIP-IN-GARCH (p, q) case. The TVZIP-INGARCH (p, q) process is discussed below, and the procedure for the TVZIP-INARCH (p) case follows in a similar manner. Expectation maximization estimation for the TVZIP-INGARCH (p, q) process. Let X 1 , X 2 ,. . .,X N be generated according to the model (2). There are two types of zeros generated by this model. They are the zeroes arising from the Poisson distribution with intensity parameter λ t and the zeroes generated by a Bernoulli process with the probability of obtaining a zero specified by the zero-inflation parameter. Therefore, a given observation can be hypothetically categorized as arising out of a Bernoulli process or as an observation from the Poisson distribution. Let us define fZ t : t 2 Zg to be a Bernoulli random variable such that Z t = 1 if X t is a generated from the Bernoulli process and Z t = 0 if it is generated by the Poisson distribution. Then, with the original parameters renamed as ϕ k , k = 1, 2,. . .,r+p+q+2. This simplified notation is used in situations where generic statements are made without reference to a specific portion of (2). Paralleling the derivations in Zhu [2], the conditional log likelihood can be written as (see S2 Appendix for details), The first derivatives of the conditional log likelihood function (7) with respect to Γ ¼ ðg 0 ; g 1 ; . . . ; g r Þ T and Θ ¼ ðy 0 ; y 1 ; . . . ; y pþq Þ T are as follows: Finally, by combining (8) and (9) the first derivative of the conditional log likelihood function with respect to F is given by: The two-step (E step and M step) Expectation Maximization algorithm is now used to estimate the parameter vector Φ ¼ ðΓ T ; Θ T Þ T . Let t t ¼ EðZ t jX t ; ΦÞ and we replace Z t byẐ t ¼ t t and define Z ¼ ðZ 1 ; Z 2 ; . . . ; Z N Þ T . Following this replacement of Z in the conditional log likelihood function, lðΦ;ẐÞ is maximized. E Step: Determine τ t using the equation Step: After Z t is replaced by its estimate, we proceed to maximize lðΦ;ẐÞ. First set dlðΦÞ d� k ¼ 0; for k = 1,2,. . .,r+p+q+2.
If,Φ, the solution to the system of equations in (10) exists, then SðΦÞ ¼ 0, where S(F) is the Fisher's score matrix, andΦ is the vector that maximizes the log likelihood, thus providing us with the estimate of the parameter vector Since a closed form solution does not exist, we require an iterative procedure to find the estimates. Let us consider the first order Taylor expansion of SðΦÞ evaluated at the valueΦ around the initial parameter values F 0 , yielding SΦ We also let the matrix of the second derivatives of the log likelihood function to be defined as From the above, we obtain the first order approximationΦ ¼ Φ 0 À H À 1 ðΦ 0 ÞSðΦ 0 Þ, and this result provides the standard Newton-Raphson algorithm. For an appropriately chosen initial valueΦ ð0Þ , the above Newton-Raphson algorithm can be used to obtain a sequence of improved estimates recursively. The improved estimates at i th iteration are updated as the initial values for the next iteration as follows: Φ ðiþ1Þ ¼Φ ðiÞ À H À 1 ðΦ ðiÞ ÞSðΦ ðiÞ Þ: This process is repeated until the differences between successive estimates are sufficiently close to zero. In our study, convergence of the EM procedure was determined by using the criterion: (2) is,

Maximum likelihood estimation for the TVZIP-INGARCH (p, q) process. The conditional likelihood function L(F) of the TVZIP-INGARCH model
The conditional log likelihood function, l(F) obtained from (11) is given by and The first derivatives of the conditional log likelihood function (12) are as follows, We can use Newton-Raphson (NR) iterative procedure to obtain the maximum likelihood estimated for the Eq (12) by setting dlðΦÞ d� k ¼ 0 for all k. With a reasonable initial starting valuê Φ ð0Þ , the i th iteration is calculated usingΦ ðiþ1Þ ¼Φ ðiÞ À H À 1 ðΦ ðiÞ ÞSðΦ ðiÞ Þ, where SΦ À � ¼ dlðΦÞ dΦ j Φ¼Φ Þ and HΦ We stop the algorithm once pre specified convergence criteria is satisfied.
Both the EM and the ML estimation procedures require careful selection of initial values of the parameter vectorΦ ð0Þ ¼ ðΘ ð0Þ ;Γ ð0Þ Þ to initiate the iterative algorithm. As described in Fokianos et al. [9], the starting valueΘ ð0Þ is computed by fitting a ARMA (p, q) model to the data. The initial value for the parameter inΓ ð0Þ is selected from a range of values on a grid that satisfies the conditions of the time-varying zero inflation function. As discussed in Weiß [31], when q>0 the parameter estimates of a INGARCH (p, q) are extremely sensitive to the choice of the initial conditional mean. In other words, inappropriate initialization ofl ð0Þ 1 may exhibit a significant effect onΘ. To address this issue, two courses of actions for computingl ð0Þ 1 are suggested. One is to treatl ð0Þ 1 as a parameter during estimation such thatl ð0Þ 1 ¼â ð0Þ 0 and the other is to use a fixed value, for examplel 1 ¼ � x. The results (see S3 Appendix) show that both initialization procedures produce relatively similar estimates, with a slightly better fit for the data observed for the initialization case withl ð0Þ 1 ¼â ð0Þ 0 . Therefore, the initial values of this study are specified according tol ð0Þ 1 ¼â ð0Þ 0 case.

Simulation study
We investigated the finite sample performance of estimators using a simulation study. The poissrnd function of MATLAB software was employed to generate the relevant data, based on recursively computed conditional intensity parameter. In order to initiate the recurve process, the intensity at time t = 0 and count data at times t�0 were set to zero (i.e., λ 0 = 0 and X l = 0 for l�0). For time periods t�1, X t was generated as follows. For each t 2 N let U t be a random variable generated from a uniform (0,1) distribution and let ω t be the zero-inflated probability at time t. Then, X t was set to zero if U t �ω t , Otherwise X t was generated from the Poisson distribution with intensity parameter λ t , where λ t was updated recursively using Eq (2). The process was repeated until the complete time series of length N was generated. Note that this is the same procedure Zhu [2] employed to generate data for his simulation study. Lengths of the time series studied were set to N = 120 and N = 360, and thousand (m = 1000) simulations runs were carried out for each parameter and sample size combination. We carried out two separate sets of simulation studies based on the two types of zero-inflation function introduced in Section 2. The profile log likelihood function given in Eq (12) was maximized using the constrained nonlinear optimization function fmincon in MATLAB. The zero-inflation probability (ω t = g(V t ,Γ)) was allowed to vary cyclically as a deterministic function of time or to be driven by an exogenous variable. Following Zhu [2], the Mean Absolute Deviation Error (MADE) was utilized as the evaluation criterion. The MADE is defined as, 1 m X m j¼1 j� j À �j where m is the number of replications and � 2 Φ ¼ ðΓ T ; Θ T Þ T is the true value while� j is the estimated value of ϕ at j th replication run. In addition, computational effort is expressed by CPU time (in seconds) and it is used as a criterion to evaluate the performance. Simulation results for Case 1 and Case 2 are reported below.
Simulation results for Case 1: Deterministic sinusoidal zero-inflation function. In this portion of the simulation study, the sinusoidal zero-inflated function ω t = g(V t ,Γ) expressed in Eq (5) was used to generate cyclically varying zero-inflation probabilities between zero and one. The following constraints are set to the parameters in the vector Γ ¼ fðA; B; CÞ T jA; B; C 2 Rg : C ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi Note that the above constraints, with δ = 0.0001, were applied in our simulation study in order to bound the zero-inflation probabilities between 0 and 1. A very small value for δ was selected to allow wider bounds for ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi , |A|, and |B|. Tables 1 through 3 provide the simulation results for the MLE estimation technique, while  Tables 4 through 6 provide simulation results for the case where estimates were obtained using the EM algorithm. The frequency of the sinusoidal wave was set at s = 12, mimicking a 12-month cycle present in monthly data. The parameter vector for the simulation study was expressed as Φ ¼ ðA; B; a 0 ; a 1 ; Note that the two estimations procedures were run on identical simulation samples for each model, parameter, and sample size combinations and hence variations due to sampling error will not be seen when comparing across estimation methods. The following tables provide summary data from the simulation runs, with the actual parameter values, estimated values, the mean absolute deviation (MADE) between them, and the CPU time taken by the estimation procedures.
The above simulation results show that the MLE and EM procedures produced almost identical means of estimates for the parameters in TVZIP-INARCH (1) and TVZIP-INARCH (2) models. For example, the means of estimates in Table 1 are almost identical to the corresponding means of estimates in Table 4. The MADE values for corresponding estimates are also almost identical across the two estimation methods. This similarity extends to means of corresponding parameter estimates across Tables 2 and 5 as well. Even in cases where the means of estimates are not identical, they are extremely close. For instance, for Model A1 with N = 120 and the true value of α 0 equals to 1.00, the mean of the MLE estimates for this parameter is 1.0473 while the mean of the EM estimates is 1.0472.
However, for the TVZIP-INGARCH (1, 1) process, there are relatively larger differences between means of estimates for the MLE and EM methods. For example, when N = 120 with true parameter values of α 0 = 1.00, α 1 = 0.20, β 1 = 0.20, the means of parameters estimates are  Table 3, Model C3), but when the sample size increases to 360, the means of the corresponding estimates changed to 0.4961 and 0.2008 respectively, which is a movement in the wrong direction. The MADE values decreased consistently with increasing sample size for both the MLE and EM estimation methods. Another relevant observation is that the simulation results for the INGARCH portion (α 1 ,β 1 ) of the model behave similar to the results obtained by Zhu [2] in the sense that the mean of the estimates are not very close to the true values even with higher sample sizes. We investigated whether this relative inaccuracy of the estimates is due to the initial values that were used for these two parameters. This was done by creating a grid of potential initial values centered on the original initial values obtained by fitting an ARAMA to the data. This approach did not improve the results and thus we reverted to using the original initial values. It is possible that many combinations ofâ andb values provide very similar likelihoods and thus the algorithm does not converge to the true values efficiently. Another important observation is that the EM algorithm took more CPU time than the MLE method across all parameter and model combinations.

Simulation study for Case 2: Zero-inflation function driven by an exogenous variable
In this part of the study, we allow the exogenous variable to generate zeros through a logistic model as described in Eq (6). The parameter vector for the simulation study under this sce-    Tables 7 through 9 provide the simulation results for   the MLE estimation techniques, while Tables 10 through 12 provide EM algorithm related simulation results.
In general, the EM algorithm takes more computational effort to estimate parameters when compared to the same model parameter estimation done by the MLE procedure. Both the MLE and EM methods produced fairly accurate estimates for the parameters across both TVZIP-INARCH (p) models. However, for the TVZIP-INGARCH (1, 1) case where the GARCH portion of the parameters (α 1 ,β 1 ) show relatively less accurate estimates, for both procedures even for the large sample size case. For instance, when N = 120 with true values of α 1 = 0.40 and β 1 = 0.30 the means of the corresponding parameter estimates are 0.3916 and 0.2439 respectively, for the MLE method, and when N = 360 these means of estimates change to 0.4794 and 0.2205 respectively (see Table 9, Model C3). A similar phenomenon is seen in the simulation results of Zhu [2]. As done in Case 1, using additional initial value combinations did not yield any improvement in this regard. In general, MADE values decrease with the increase in sample size.

Extension to the softplus INGARCH model
Herein we introduce a generalization of the softplus INGARCH model proposed by Weiß et al. [12] by incorporating zero-inflation to the underlying count-data distribution and allowing a set of exogenous variables to act as inputs to the model that defines the conditional mean of the process. In the following, we assume that the underlying distribution is zero-inflated versions of the Poisson or the negative binomial distributions or the generalized Poisson. The zero-inflated softplus INGARCHX model we developed is defined as follows. Let fX t : t 2 Zg is the count process, and μ t defines the conditional mean of X t given the sigma-fields F X tÀ 1 and F V t such that m t ¼ EðX t jF X tÀ 1 ; F V t ; ΘÞ. Here, F X tÀ 1 is the sigma-field generated by fX l ; l l ; l � t À 1g and F V t is the sigma-field generated by the set of exogenous covariates {V 1 ,. . .,V t } with V j ¼ ðV j1 ; V j2 ; . . .; V jr Þ T for j = 1,. . .,t. Note that Θ is a set of parameters defined later in the section. Let qð:jl * t ; ΩÞ be a probability mass function of a discrete random variable U t U t , with conditional mean l * t , for t 2 Z. We define λ t such that l t ¼ l * t if the distribution of U t is Poisson or negative binomial, and l t ¼ ð1 À φÞl * t if the distribution is generalized Poisson, where φ is the dispersion parameter. This defines the underlying count data distribution without zero inflation. Now let ω be the zero-inflation probability associated with X t . Then the conditional probability mass function of X t jF X tÀ 1 ; F V t ; Θ is given by,

PLOS ONE
for t 2 Z. If qð:jl * t ; ΩÞ is representing a Poisson, negative binomial Type 1 or negative binomial Type 2 distribution, then the g(λ t ) = λ t . If the underlying distribution is generalized Poisson, then gðl t Þ ¼ l t 1À φ : If the proposed softplus INGARCHX version, we let λ t to be defined as with the softplus link function given by s c z ð Þ ¼ cln 1 þ exp z c À � À � with c>0. The set of parameters in the model can now be defined as Θ ¼ ðo; a 0 ; a 1 ; . . . ; a p ; b 1 ; . . . ; b q ; g 11 ; . . .; g tr ; c; ΩÞ:

Real data examples
In this section the proposed TVZIP-INGARCH (p, q) model is applied to two real-world datasets. The time-varying component given by Eq (5) can be formulated in many ways. In example one, we selected two scenarios for this component, and labeled them Sc1 and Sc2. In example two, the time-varying component was labeled as obeying Scenario 4 (Sc4). Details  [21], were also fitted to the data. In addition, the following formulations were also included in the comparison. The log-linear INGARCH model of Fokianos and Tjøstheim [10] was modified to include zero-inflation. We also modified the log-linear INGARCHX model of Chen and Lee [11] to include zero-inflation, but it should be noted that with the inclusion of zero-inflation the above two models are nested within the non-time-varying version of the model of Xu et al. [30]. We extended the softplus INGARCH model introduced by Weiß et al. [12], as described in an earlier section, by introducing zero inflation and an exogenous variable.  [2]. Note that the compound Poisson model of Gonçalves et al. [21] was fitted assuming the zero-inflated geometric Poisson and the zero-inflated Neyman Type A distributions. The above combinations of model categories and distributions yield a total of twentythree count data time series formulations. These twenty-three formulations are listed in Table 1 in S4 Appendix. These formulations were then fitted using three model orders, namely M1 with p = 1,q = 0; M2 with p = 2,q = 0; and M3 with p = 1,q = 1. This gives rise to sixty-nine total models. Not all of these sixty-nine models were fitted to data in each of the two examples shown below. Example 1 does not contain data on an exogenous variable and hence the log-linear INGARCHX and the softplus INGARCHX model categories were not used in the first example. Models belonging to the above two categories, however, were fitted to the data in Example 2 in place of the log-linear INGARCH and the softplus INGARCH model categories. Excluding versions of the proposed model, forty-five models with constant zero-inflation probability were fitted to each of the example data sets, with eighteen models common to both sets (twelve models from C1 category and six models from C2 category).
Since we are comparing the proposed time-varying zero inflation model against models with constant zero-inflation probability, an argument can be made that any superior performance by the proposed model can be due to the presence of a change point. Thus, fitting a suitably selected constant zero-inflation model separately to the two sub-series before and after a detected change point may negate any such superiority. To address this potential criticism, the following approach was taken with respect to the INGARCH models (C1 category) with zeroinflated versions of Poisson, negative binomial Type 1, negative binomial Type 2, and the generalized Poisson distributions under model orders M1, M2, and M3. First, we fitted the above set of constant zero-inflation models to the example data. Then the model with the smallest AIC value from among the zero-inflated INGARCH models within each model order M1, M2, M3 was selected. We repeated the same procedure using the Bayesian Information Criteria (BIC) and found that the models selected using the AIC values also had the lowest BIC values. Following that, the CUSUM test proposed by Lee et al. [27] was applied using the residuals from the selected model. If a change point was detected, the model under consideration was fitted independently to the two subseries before and after the change point. It is noteworthy that a change point was detected every time this test was carried out. The AIC and BIC values were computed for all models, including in cases where the model was fitted to two sub-series. Note that the theoretical properties of the CUSUM test proposed by Lee et al. [27] are based on ten assumptions which holds for the zero-inflated versions of Poisson, negative binomial, and generalized Poisson, according to Lee et al. [27] and Lee and Lee [32]. We did not, however, verify whether models in other categories would also satisfy all these assumptions. Thus, we carried out this procedure for all models falling into the zero-inflated INGARCH (C1) category but did not do so for other model categories.
Following the above process, we identified models with the lowest AIC and BIC values within each model category (C1, C2, C3, C4, C5, C6) by model order (M1, M2, M3) combination. Models thus identified can be seen in cells highlighted in light gray in Table 1B of S5 Appendix (for Example 1) and Table 1B of S6 Appendix (for Example 2). Note that Table 1A in each Appendix provides the AIC and BIC values for versions of our proposed model. Within each model category (but across model orders), we identified the model with the lowest AIC value and the model with the lowest BIC value among models highlighted in light gray. Then from among the above models, we also identified the best model(s) across all model categories by using AIC and BIC values. After the best model or models across all model categories were selected, standardized Pearson residuals were calculated, and used to check model adequacy.
The first example is based on the Influenza A associated pediatric deaths data set downloaded from the Centre for Disease Control and Prevention web page. This data was modelled by using the TVZIP-INGARCH (p, q) formulation with a sinusoidal zero-inflation function. In the second example, we used the Pediatric mortality caused by Influenza B data set downloaded from the same web page to demonstrate the performance of TVZIP-INGARCH (p, q) process where the zero-inflation is driven by an exogenous variable. The data sets used in this study are available in the supplement labeled S1 Data.

Real data example-Use of a deterministic sinusoidal zero-inflation function
The TVZIP-INARCH (1), TVZIP-INARCH (2), and TVZIP-INGARCH (1, 1) models with sinusoidal zero-inflation were fitted to the Influenza A associated pediatric mortality data set. Performance of the proposed models was compared to those of the constant zero-inflated count data time series models mentioned above.
The data set provides weekly count data of U.S. pediatric deaths caused by type A influenza virus over the time period from week 40 of 2015 to week 43 of year 2018. This consists of 160 weekly observations of pediatric death counts. The data were extracted from weekly U.S. Influenza Surveillance report, which was published by Centre for Disease Control and Prevention (CDC) [33]. Summary statistics of the data showed a mean of 1.5063 and a variance of 6.352, suggestive of overdispersion. Fig 1 illustrates the frequency of each pediatric mortality case caused by influenza A virus using a bar chart. Observe that there are 86 zeros, which comprises 53.75% of the total time points. Fig 2 illustrates the original time series of the counts followed by the plots of its autocorrelation function (ACF) and the partial autocorrelation function (PACF), respectively. The bar chart shows the excess number of zeros in the data, and the time series plot demonstrates prolonged periods of zero counts, which supports the use of zero-inflated Poisson time series models to analyze this data. Furthermore, we can observe an annual seasonality in the portions of the time series with excessive zero counts.
To understand the zero-inflation behaviour of this data set, we aggregated weekly data in to its corresponding calendar month and constructed the total monthly zero mortality counts. These counts were then converted into a monthly proportion by dividing each monthly total

PLOS ONE
zeros by the maximum zero count over the observed period to scale values at or below one. According to the Fig 3, the plot of the monthly proportion of zero counts exhibits an approximate sinosidual behaviour throughout the observed time span. Thus, the general sinosidual function mentioned in Eq (5) can be used to model the zero-inflation behaviour of this data. More details of this modelling will be described later in this section.
As mentioned earlier, two scenarios were considered for the time-varying zero-inflation component of the proposed model. In the models listed under Scenario 1 (Sc1) we assumed a piecewise constant zero-inflation fucntion. The monthly zero-inflation was allowed to vary according to a sinusoidal zero-inflation function with a 12-month period, but the weeks within a given month were assigned the same zero-inflated probability value associated with that month. The models under Scenario 2 (Sc2) used a sinusoidal function describing a time-

PLOS ONE
varying zero-inflation probability that changes on a weekly basis. Note that all other model categories we consider fall under Scenario 3 (Sc3) where the zero-inflation probability remains constant over time. In addition to the proposed model under scenarios Sc1, and Sc2, we fitted models in categories C1, C2, C3, and C5. All models were fitted under the model orders M1, M2, M3 with various underlying zero-inflated distributions described at the beginning of the real data example section. The EM algorithem was used to estimate the model parameters and the results for Sc1 and Sc2 scenarios are reported in the Table 13. The EM algorithm was chosen over the MLE method to be in line with the method used in Zhu [2] for its real data example.  Table 13 shows that the models which fall under the Sc1 exhibited lower AIC and BIC values compared to their counterparts under Sc2. We could conclude that the use of a piecewise constant time-varying zero-inflation function improved the fit compared to version where the time-varying zero-inflation probability changes on a weekly basis.

PLOS ONE
As described ealier, both the AIC and BIC values were used to identify the best models within each model category by model order combination for the rest of the constant zeroinflation models. The models thus identified are shown in light gray highlighted cell in Tables 1A and 1B in S5 Appendix. Models with the lowest AIC or BIC values within each model category are identifies in boldface type in Tables 1A and 1B in S5 Appendix. If the model has the lowest information criteria values with respect to both AIC and BIC, then it is identified in italic boldface font. These models with the lowest AIC or BIC within each model category are also shown listed in Table 14, togther with the Sc1 and Sc2 versions of the proposed model under model order M2, which was selected as the best among all three orders.
The performance of the two versions of the TVZIP-INGARCH (p, q) models with the lowest AIC and BIC values are compared with those of other INGARCH time series models that assume constant zero-inflation in Table 14. Among all the models, the TVZIP-INARCH (2) under Sc1 zero-inflation assumption has lower AIC and BIC values. Note that the second best model with respect to AIC and BIC criteria is also a version of the proposed TVZIP-IN-GARCH model.
The TVZIP-INARCH (2) model, with the piecewise zero-inflation formulation, which has the lowest AIC and BIC values amoung all models considered, was selected for further investigation. The standardized Pearson residuals were calculated for the selected model and used for diagnostic tests. The ACF plot and the histogram, both based on the standardized Pearson residuals are presented in  The standardized Pearson residuals display a mean (0.0017), which is close to zero, and a variance (0.8108), which is close to one. The ACF plot indicates that residuals of the fitted model meet the assumption of no autocorrelation at 5% significance level. Additionally, the mean of the fitted model (1.4608) is close to the observed mean of the count data series (1.5063). Based on results from both model selection criteria, model fit, and residual diagnostic results, we conclude that the TVZIP-INARCH (2) model with Sc1 formulation for zero-inflation probability provided the best fit to the data. In this model, we assumed that weeks within any given month has a constant zero-inflation, yet monthly zero-inflation varies cyclically.

Real data example-Zero-inflation function is driven by exogenous variable
In this subsection we examine the performance of the TVZIP-INGARCH (p, q) models when the zero-inflation function is driven by an exogenous variable. Influenza B associated pediatric mortality data set was used, and the TVZIP-INARCH (1), TVZIP-INARCH (2), and TVZI-P-INGARCH (1, 1) models were fitted to the data. Results were compared with forty-five constant zero-inflated INGARCH (p, q) models. We selected the weekly average of nationwide low temperatures as the exogenous variable that drives the zero-inflation probability. This is motivated by one of the reasons generally accepted as a cause for the easy transmission of influenza during winter months, namely the cold temperatures driving people indoors. Influenza B associated pediatric mortality data set was accessed from the weekly U.S. Influenza Surveillance Report [33], and the temperature data were extracted from the United States National Oceanic and Atmospheric Administration (NOAA) weather prediction center [34]. Both data sets spanned the time period from week 40 in year 2014 to week 39 in year 2018. The influenza B data set contained 209 weekly observations of pediatric death counts. Summary statistics of pediatric mortality cases due to influenza B showed a mean of 0.8517 and a variance of 2.4538. Since the empirical variance was higher than the empirical mean, the data exhibit overdispersion. Fig 5 illustrates the frequency of each pediatric mortality count caused by influenza B virus type using a bar chart. The bar chart shows that there are 128 zeros, which comprises 61.24% of total of the time points. This suggests that the pediatric mortality data are zero-inflated.
The time series plot, ACF, and PACF plots are given in Fig 6. Based on time series plot, we can see that there is an annual seasonality exhibited in this data set. Moreover, it shows that

PLOS ONE
there were periods with clusters of zeros that occur periodically. Therefore, as discussed in the above section, we suggested TVZIP-INGARCH (p, q) to model the count data series.
In this example, the time-varying zero-inflation was modelled by considering an exogenous time series, namely the weekly average of nationwide low temperature. Comparison of the two time series plots is given in Fig 7. Fig 7 shows that periods with higher values of low temperature coincide with periods of zero pediatric mortality caused by Influenza B. Hence, it is suggestive that periods with high zero counts (low pediatric mortality) are related to periods with higher values of low temperature. Thus, we used the weekly average low temperature data to model the time-varying zero-

PLOS ONE
inflated component in the pediatric mortality data set. In this example, time-varying zeroinflation process was modelled using the formulation from Eq (6) with low temperature as the input variable to the logistic model. We label this zero-inflation scenario as Sc4. Results from this the TVZIP-GARCH formulation under model orders M1, M2, and M3, with varying underlying distributions are given in Table 15. Note that we used the EM algorithm as was done in Example 1.
Results in Table 15 show that the TVZIP-INARCH (2) model fitted under Scenario 4 (Sc4) exhibited the lowest AIC and BIC values compared to the other models. Forty five models were fitted under the constant zero-inflation scenario (Sc3). Table 1B in S6 Appendix provides list of these models, which fall into categories C1, C2, C4, and C6 with model orders M1, M2, and M3. The zero-inflated versions of Poisson, negative binomial Type 1, negative binomial Type 2, and generalized Poisson distributions we considered as the underlying distributions associated witn the counts in categories C1, C4, and C6. In C2, we assumed the zero-inflated geometric Poisson and zero-inflated Neyman Type A as the underlying distributions. As was done in Example 1, a test for a change point was conducted if a model had the lowest AIC value under each model category C1 by model order combination. Note that at this initial stage, models with the lowest AIC value under the C1 category had the lowest BIC value as well. If such a point was detected, the model under consideration was independently fitted to sub-series before and after the change point. Comparison of the versions of the proposed model with the forty-five existing models, as well as for models fitted to subseries, are repoted in Table 1A and 1B in S6 Appendix. The models with the lowest AIC or BIC values in each model category by model order combination are highlighted in light gray. Models with the lowest AIC or BIC values within each model category are identifies in boldface type in Tables 1A and 1B in S6 Appendix. If the model has the lowest information criteria values with respoect to both AIC and BIC, then it is identified in italic boldface font.
The results for TVZIP-INGARCH (p, q) models under Scenario 4 (Sc4) together with the contant zero-inflation models under Scenario 3 (Sc3) that had the lowest AIC or BIC values within each model category are shown in Table 16.
Based on the BIC criterion, TVZIP-INARCH (2) model provided the best fit to the data, albeit by a marginally lower BIC value compared to the next best model. However, a different result is obtained when considering the AIC values. AIC values indicate that the ZINB2 softplus INARCHX (2) model, which assumes a constant zero-inflation, exhibits a marginally better fit to the data. This is the new model we introduced as a generalization of the softplus INGARCH formulation of Weiß et al. [12] by including zero-inflation and an exogenous variable to the original version. The inclusion of an exogenous variable seems to lead to an enhanced performance with respect to AIC. This advantage vanishes slightly when BIC is used as a model selection criterion. Since both models perform well, model adequacy checks were

PLOS ONE
carried out for the TVZIP-INARCH (2) and the ZINB2 softplus INARCHX (2) models. The ACF plots, and the histograms of standardized Pearson residuals for the fitted models are presented in Figs 8 and 9 respectively. The ACF plots indicate that the standardized Pearson residuals of both TVZIP-INARCH (2) and the ZINB2 softplus INARCHX (2) models do not exhibit significant autocorrelations. However, the mean (-0.0054) and the variance (1.0006) of the residuals of the TVZIP-I-NARCH (2) model are closer to 0 and 1 respectively. The residuals of the ZINB2 softplus INARCHX (2) model have a mean (-0.0389) close to 0, but variance (0.8557) is further away from 1 compared to the variance of the TVZIP-INARCH (2) Model. While both models satisfied the residual checks for model adequacy, the TVZIP-INARCH (2) may be considered to perform marginally better due to its slightly lower BIC value and the residual mean and variance being closer to the ideal values of 0 and 1 respectively. In addition to that, the mean of the fitted values of TVZIP-INARCH (2) model (0.8421) is closer to the sample mean (0.8517) than that of the ZINB2 softplus INARCHX (2) model (0.9264). Thus, we can conclude that while both these models performed well, our proposed TVZIP-INGARCH model has a tenuous edge over the Softplus INGARCHX model. In a practical sense, however, one can argue that both models perform equally well.
One reason for the good performance of ZINB2 softplus INARCHX model may be due to the nature of the softplus function that can bring the value of the conditional mean close to zero for some values of the exogenous variable. Thus, while the zero-inflated component was set up to be non-time-varying, the influence of the exogenous variable can bring the conditional mean of the process very close to zero over some periods. In other words, softplus link function allows the exogenous variable to significantly increase the probability of obtaining zero counts in a time-varying fashion.

PLOS ONE
The main aim of this example is to illustrate the utility of the proposed formulation in modelling a real-life situation. Based on the results from this example, TVZIP-INARCH (2) formulation with zero-inflation modeled using an exogenous variable and the ZINB2 softplus INARCHX (2) model appear as good candidates to be included in the tool set a practitioner may consider when modeling count data time series with time-varying zero inflation, in situations where an appropriate exogenous variable is available. Equally important is the observation that no other model outperformed the proposed model with respect to any criteria based on residual analysis.

Conclusions
A time-varying zero-inflated Poisson integer GARCH (TVZIP-INGARCH) model was proposed to accommodate situations where zero-inflation is driven by either a deterministic function of time or a set of exogenous variables. Monte-Carlo simulation results indicate that the Expected Maximization (EM) and Maximum Likelihood Estimation (MLE) methods produce very similar results with respect to parameter estimates. It is observed that both EM and MLE techniques estimated the model parameters with good accuracy when the underlying model has a purely ARCH component. When the model has a GARCH component, the GARCH parameters are estimated with lesser accuracy, which is a phenomenon also seen in the study of the existing ZIP-INGARCH model proposed by Zhu [2]. When tested on two real-life data sets, the TVZIP-INGARCH models performed better than the other constant zero-inflated INGARCH formulations with one possible exception. Even in this case the competing model does not outperform the proposed model. These results illustrate the usefulness of the proposed model as one of the tools the practitioner can utilize for modeling empirical count data time series with time-varying zero-inflation. In addition, the flexibility of modelling zero-